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Abstract 

This  report  describes  some  of  the  problems,  achievements,  and  directions  of  investigation  of  three  ar¬ 
eas  of  research.  The  first  is  the  development  of  combinatorial  optimization  techniques  to  solve  the  central 
problem  of  multisensor  and  multitarget  tracking,  i.e.,  the  data  association  problem  of  partitioning  obser¬ 
vations  into  tracks  and  false  alarms.  The  problem  formulation,  algorithm  design,  and  real-time  solution 
techniques  involve  from  probability  and  information  theory,  system  identification,  filtering,  control  systems, 
combinatorial  optimization,  and  advanced  computer  architectures,  including  massively  parallel  computers. 
The  data  association  problem  for  general  multitarget/multisensor  tracking  problems  is  posed  as  a  class  of 
multidimensional  assignment  problems.  The  algorithms  under  development  are  based  on  a  recursive  La- 
grangian  relaxation  scheme,  construct  near-optimal  solutions  in  real-time,  and  use  a  variety  of  techniques 
ranging  from  two  dimensional  assignment  algorithms,  a  conjugate  subgradient  method  for  the  nonsmooth 
optimization,  graph  theoretic  properties  for  problem  decomposition,  and  a  branch  and  bound  technique  for 
small  solution  components.  A  model  problem  is  presented  to  demonstrate  the  efficiency  and  robustness 
of  the  current  algorithms.  The  second  part  centers  on  investigation  of  various  numerical  methods  for  the 
solution  of  nonlinear  optimal  control  problems.  The  analysis  of  convergence  in  infinite  dimensional  spaces, 
discretizations,  and  numerical  implementations  are  in  progress  for  Newton,  penalty,  augmented  Lagrangian, 
and  interior/exterior  point  methods.  The  final  part  is  the  investigation  of  parametric  constrained  optimiza¬ 
tion  problems  using  numerical  bifurcation  and  continuation  methods  with  applications  to  design  optimization 
and  control  systems. 
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I.  Introduction 


This  report  describes  some  of  the  problems,  achievements,  and  directions  of  investigation  of  three  areas 
of  research  as  described  below.  Over  the  past  two  years  twenty  three  presentations  have  been  given  and 
eighteen  papers  have  been  published  or  submitted  for  publication,  2  MS  and  4  PhD  students  have  graduated, 
and  A.  B.  Poore  has  been  appointed  associate  editor  of  Computational  Optimization  and  Applications.  To 
reorient  this  program  to  the  area  of  multisensor  data  fusion,  Professor  Poore  spent  the  summer  at  Rome 
Labs  as  part  of  the  AFOSll  Summer  Faculty  Research  Program  as  described  in  Section  II  The  technical 
information  for  the  last  two  years  is  given  in  Section  VII.  The  three  areas  of  research  are  briefly  explained 
in  the  remainder  of  this  introduction. 

The  first  part  of  this  research  program  centers  on  the  development  of  combinatorial  optimization  tech¬ 
niques  to  solve  the  central  problem  of  mulitsensor  data  fusion  and  multitarget  tracking,  i.e.,  the  data  associ¬ 
ation  problem  of  partitioning  observations  into  tracks  and  false  alarms.  The  problem  formulation,  algorithm 
design,  and  real-time  solution  involve  techniques  from  probability  and  information  theory,  system  identifica¬ 
tion,  filtering,  control  systems,  combinatorial  optimization,  and  advanced  computer  architectures,  including 
massively  parallel  computers.  The  data  association  problem  for  general  multitarget  tracking  problems  is 
formulated  as  a  class  of  multidimensional  assignment  problems.  The  algorithms  under  development  are 
based  on  a  recursive  Lagrangian  relaxation  scheme,  construct  high  quality  suboptimal  solutions  in  real-time, 
and  use  a  variety  of  techniques  ranging  from  two  dimensional  assignment  algorithms,  a  conjugate  subgra¬ 
dient  method  for  the  nonsmooth  optimization,  graph  theoretic  properties  for  problem  decomposition,  and 
a  branch  and  bound  technique  for  small  solution  components.  This  problem  of  partitioning  multiple  data 
sets  at  some  cost  or  to  some  benefit  is  also  the  central  problem  in  perceptual  grouping  in  psychology  and 
stereo  correspondence  in  both  biological  and  computer  vision  [9].  Thus  the  applications  potentially  extend 
far  beyond  the  current  applications.  The  current  status  and  results  of  this  research  effort  are  described  in 
Section  III. 

The  second  part  centers  on  the  investigation  of  various  numerical  methods  for  the  solution  of  nonlinear 
optimal  control  problems.  The  analysis  of  convergence  in  infinite  dimensional  spaces,  discretizations,  and 
numerical  implementations  are  in  progress  for  .Newton,  penalty,  augmented  Lagrangian,  and  interior  point 
methods.  A  longer  term  goal  is  the  investigation  of  parametric  problems  in  nonlinear  control  systems 
including  but  not  limited  to  the  nonlinear  optimal  control  problem.  Some  of  the  initial  results  in  this 
direction  are  described  in  Section  IV. 

The  final  part  of  this  research  program  is  the  investigation  of  parametric  nonlinear  programming  prob¬ 
lems  using  numerical  bifurcation  and  continuation  methods  with  applications  to  design  optimization  and 
parametric  control  systems,  and  represents  a  potential  for  a  real  extension  of  our  understanding  of  basic 
phenomena,  global  sensitivity,  robustness,  and  multiplicity  of  solutions  in  much  the  same  way  that  these 
theoretical  and  numerical  techniques  have  helped  the  understanding  of  dynamical  systems  and  nonlinear 
equations.  Thus  the  objective  in  this  aspect  of  the  research  program  is  to  develop  the  analytical  and  numeri¬ 
cal  techniques  to  map  out  regions  of  qualitatively  different  behavior  and  to  locate  the  “stability”  boundaries 
of  these  regions  in  parameter  space.  The  latter  is  important  because  drastic  changes  in  the  optimum  occur 
in  the  presence  of  singularities  which,  in  turn  defu„-  ibosc  “ct-iKiiity”  KrmnfWipc  ^uch  Knowledge  allows  for 
the  uncertainty  in  system  and  model  parameters  and  yields  information  about  the  expected  behavior  when 
control  parameters  are  varied  to  enhance  the  performance  of  the  system  under  consideration.  In  addition  to 
providing  a  global-like  sensitivity  analysis,  these  methods  are  quite  efficient  in  computing  multiple  optima. 
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Several  model  problems  taken  from  the  very  active  area  of  design  optimization  are  being  investigated  to  test 
and  illustrate  the  value  and  applicability  of  these  continuation  and  bifurcation  methods,  as  well  as  to  provide 
motivation  and  focus  for  lurcher  development.  Preliminary  theory  and  numerical  implementation  have  been 
completed  as  described  in  detail  in  Section  V. 

II.  AFOSR  Summer  Faculty  Research  Program  at  Rome  Labs,  Griffiss,  AFB 

During  the  past  five  years  we  have  worked  on  multitarget  tracking  arid  less  on  multisensor  data  fusion. 
Thus  to  reorient  our  research  program  more  in  line  w'ith  the  needs  of  the  Department  of  the  Air  Force  in 
multisenor  data  fusion,  1  spent  the  summer  at  Rome  Labs  as  part  of  the  AFOSR  summer  faculty  research 
program.  My  goals  were:  (1)  to  understand  the  multisensor  data  fusion  and  multitarget  tracking  problem 
from  the  viewpoints  of  the  researchers  at  Rome  Labs,  (2)  to  survey  the  existing  methods  and  approaches 
to  multisensor  data  fusion  and  multitarget  tracking,  (3)  to  work  on  modeling  multisensor  data  fusion  and 
multitarget  tracking  problems,  and  (4)  to  develop  a  working  relation  with  researchers  at  Rome  Labs  with 
a  long  term  goal  of  cooperative  research.  VV'orking  with  Martin  Liggins  and  Vincent  Vannicola  at  Rome 
Labs,  we:  (1)  achieved  a  new  and  broadened  view  of  the  needs  in  multisensor  and  multitarget  tracking,  (2) 
developed  a  unifying  approach  to  muitisensor  data  fusion  and  multitarget  tracking  based  on  multidimensional 
assignment  problems,  (3)  worked  on  a  paper  that  is  to  appear  [38]  and  will  be  presented  at  the  1993  SP1E 
meeting  in  Orlando  in  April,  1993,  (4)  arrived  at  a  more  unified  understanding  of  the  elements  needed  to 
smoothly  transition  6.1  research  into  6.2  and  6.3A  areas,  and  (5)  examined  the  aspects  of  involving  Rome 
Laboratory  personnel  in  developing  an  in-house  6.1  research  base. 

III.  Optimization  Problems  in  Multitarget/Multisensor  Tracking 

This  section  describes  some  of  the  optimization  problems  in  multisensor/multitarget  tracking.  Section 
A  presents  the  problem  area  overview,  Section  B  summarizes  the  achievements,  Section  C  presents  a  case 
study,  and  Section  D  gives  an  overview  of  the  algorithms. 

III.  A.  Problem  Statement.  The  central  problem  in  any  multitarget/multisensor  surveillance  system 
is  the  data  association  problem  of  partitioning  the  observations  into  tracks  and  false  alarms.  Current  meth¬ 
ods  [7]  for  mutitarget  tracking  generally  fall  into  two  categories:  sequential  and  deferred  logic.  Methods 
for  the  former  include  nearest  neighbor,  one-to-one  or  few-to-one  assignments,  and  all-to-one  assignments 
as  in  the  joint  probabilistic  data  association  (JPDA)  [3],  For  track  maintenance,  the  nearest  neighbor 
method  is  valid  in  the  absence  of  clutter  when  there  is  no  track  contention,  i.e.,  when  there  is  no  chance  of 
misassociation.  Problems  involving  one-to-one  or  few-to-one  assignments  are  generally  formulated  as  (two 
dimensional)  assignment  or  multi-iissignment  problems  for  which  there  are  some  excellent  optimal  algo¬ 
rithms.  This  methodology  is  real-time  but  can  result  in  a  large  number  of  partial  and  incorrect  assignments, 
particularly  in  dense  or  high  contention  scenarios,  and  thus  incorrect  track  identification.  The  difficulty  is 
that  decisions,  once  made,  are  irrevocable,  so  that  there  is  no  mechanism  to  correct  misassociations.  The 
use  of  all  observations  in  a  scan  (e  g.,  JPDA)  to  update  a  track  moderates  the  misassociation  problem  and 
has  been  successful  for  tracking  a  few  targets  in  dense  clutter. 

Deferred  logic  techniques  consider  several  data  sets  or  scans  of  data  all  at  once  in  making  data  associ¬ 
ation  decisions  At  one  extreme  is  batch  processing  m  which  all  observations  (from  all  time)  are  processed 
together,  but  this  is  computationally  too  intensive  for  real-time  applications.  The  other  extreme  is  sequen¬ 
tial  processing.  Deferred  logic  methods  between  these  two  extremes  are  of  primary  interest  in  this  work. 
The  principal  deferred  logic  method  used  to  track  large  numbers  of  targets  in  low  to  moderate  clutter  is 
called  multiple  hypothesis  tracking  (MHT)  in  which  one  builds  a  tree  of  possibilities,  assigns  a  likelihood 
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score  based  on  Bayesian  estimation,  develops  an  intricate  pruning  logic,  and  then  solves  the  data  association 
problem  by  explicit  enumeration  schemes. 

Another  important  aspect  in  surveillance  systems  is  the  growing  use  of  multisenor  data  fusion  in  which 
one  associates  reports  from  multiple  sensors  together.  Once  matched,  this  more  varied  information  has  the 
potential  to  greatly  enhance  target  identification  and  state  estimation  (7j.  The  central  problem  is  that  of 
data  association  and  the  principal  method  employed  is  the  deferred  logic  technique  of  multiple  hypothesis 
tracking  [1,2, 7, 8].  This  problem  of  partitioning  multiple  data  sets  at  some  cost  or  to  some  benefit  is  also 
the  centra!  problem  in  perceptual  grouping  in  psychology  and  stereo  correspondence  in  both  biological  and 
computer  vision  [9j. 

Thus  data  association  and  the  more  general  the  problem  of  partitioning  multiple  data  sets  to  some 
benefit  is  a  fundamentally  important  combinatorial  optimization  problem.  These  problems  generally  have 
the  following  characteristics:  the  problems  are  large  scale;  the  objective  function  is  noisy  due  to  plant 
noise,  errors  in  the  sensor  measurements,  and  modeling  uncertainty;  and  they  NP-hard  [14]  but  must  be 
“solved”  real-time.  Consider  the  methods  currently  used  to  solve  these  data  association  problems:  explicit 
enumeration  and  greedy  algorithms.  The  former  is  inevitably  faulty  in  dense  scenarios  since  the  time  required 
to  solve  the  problem  optimally  can  grow  exponentially  in  the  size  of  the  problem.  The  latter  cannot  produce 
near-optimal  solutions  with  any  robustness. 

Hence  the  challenge  is  the  design  of  algorithms  that  solve  these  problems  to  the  noise  level  of  the  problem 
in  real-time.  This  is  precisely  the  achievement  of  this  research  program.  All  of  these  problems  have  been 
formulated  as  multidimensional  assignment  problems,  and  a  new  class  of  algorithms  based  on  Lagrangian 
relaxation  has  been  developed  to  construct  near-optimal  solutions  in  real-time.  The  potential  implications 
are  the  orders  of  magnitude  improvement  in  the  speed  of  existing  MHT  algorithms  and  the  extension  of  the 
problems  solvable  by  MHT. 

Finally,  it  is  important  to  note  that  the  use  of  combinatorial  optimization  in  multitarget  tracking  is  not 
new  and  dates  back  to  the  mid-sixties  and  the  pioneering  work  of  Sittler [46],  who  used  maximum  likelihood 
estimation  to  evaluate  all  possible  track  updates  and  employed  track  splitting  (several  hypotheses  were 
maintained  for  each  track)  and  pruning  (when  their  probabilities  fell  below  certain  threshold).  Maximum 
likelihood  estimation  was  further  investigated  by  Stein  and  Blackman  [6],  who  developed  a  comprehensive 
probability  for  track  initiation,  track  length  expectancy,  missed  detections  and  false  alarms.  Morefield  [25] 
pioneered  the  use  of  the  integer  programming  to  solve  a  set  packing  problem  arising  from  a  data  association 
problem.  Multiple  hypothesis  tracking  has  been  popularized  by  the  fundamental  work  of  Reid  [45].  The 
work  here  now  makes  these  approaches  practical. 

III.B.  Achievements.  Our  activities  and  achievements  in  1991  and  1992  have  been  many  and  varied 
and  can  be  briefly  summarized  as  follows: 

•  As  a  participant  in  the  AFOSR  Summer  Faculty  Research  Program,  Professor  Poore  spent  the  summer 
of  1992  at  Rome  Labs  at  Grifliss  AFB.  A  major  accomplishment  was  the  demonstration  the  existing 
deferred  logic  association  techniques  such  as  multiple  hypothesis  tracking  that  is  the  technique  used 
in  multisensor  data  fusion  and  multitarget  tracking  can  be  replaced  by  multidimensional  assignment 
problems.  This  work  is  the  suuject  of  two  forthcoming  papers  by  Drs.  A.  B.  Poore  and  N.  Rijavec  of 
CSU  and  Dr.  V.  C.  Vannicola  and  M.  Liggins  [37,38]  of  Rome  Labs. 

•  One  class  of  algorithms  for  the  construction  of  real-time  solutions  of  the  mulitsensor/muititarget  tracking 
problems  has  been  developed.  The  basic  scheme  currently  [29,40,41]  employs  preprocessing  in  the  form 


of  gating  and  clustering.  Then  the  sparse  decomposed  problems  are  solved  by  a  recursive  Lagrangian 
relaxation  scheme.  A  AT-dimensional  assignment  problem  is  relaxed  to  a  {K  —  l)-dimensional  one  by 
incorporating  one  set  of  constraints  into  the  objective  function  using  a  Lagrangian  relaxation  of  this 
set.  Given  a  solution  of  the  { K  -  l)-dimensional  problem,  a  feasible  solution  of  the  /(-dimensional 
problem  is  then  reconstructed.  The  (K  -  l)-dimensional  problem  is  solved  in  a  similar  manner  and  the 
process  is  repeated  until  one  reaches  the  twodimensional  problem  which  is  solved  exactly.  The  duality 
gap  in  this  process  is  generally  small  and  one  obtains  good  bounds  on  'he  optimal  solution.  The  full 
technical  description  can  be  found  in  the  forth  coming  paper  of  Poore  and  Rijavec  [29,40],  Other  classes 
of  algorithms  are  underdevelopment. 

•  Algorithms  for  an  initial  tracking  system  have  been  developed.  This  includes  problem  formulation 
for  track  initiation  and  track  maintenance.  Algorithms  for  system  optimization  and  estimation  (least 
squares,  Kalman  filtering,  nonlinear  estimation  techniques),  model  simulation,  solution  quality  mea¬ 
surements  have  been  developed  and  implemented.  This  provides  the  basis  for  checking  the  quality  for 
different  algorithms  under  development. 

•  Extensive  simulations  have  been  performed  to  demonstrate  the  robustness  and  speed  of  the  assignment 
solvers.  These  simulations  are  the  subject  of  previous  and  forthcoming  publications  [29-34,39,41]. 

•  Twenty  three  presentations  have  been  given  and  eighteen  papers  have  been  submitted  or  published. 
III.C.  A  Case  Study 

In  this  section,  a  model  problem  is  formulated  and  solved  to  demonstrate  the  overall  performance  of  the 
algorithms  involved. 

III.C. 1.  The  Model  Problem.  One  generally  assumes  that  each  target  is,  except  in  a  maneuver, 

modeled  by  a  state  state-space  system  of  the  form 

z(k+l)  =Fk(x(k))  +  Gk(x(k))xv(k) 


z(k)  =Hk{x{k))  +  v(k) 


where  x(k)  is  a  vector  of  n  state  variable,  tv  and  v  are  independent  white  noise  sequences  of  normal  random 
variables,  z{k)  represents  the  measurement  at  time  k  associated  with  this  particular  target  and  Hk(x(k)) 
relates  the  state  x(k)  to  the  measurement  x(k).  In  this  case  study  the  targets  are  assumed  to  travel  in  two 
dimensional  space  according  to  the  constant  acceleration  model 


,  .  t* 

x{t,Q )  -  x0  +  lvx  +  —ax 

L? 

0  =  Vo  +  tvy  +  —ay 


(C.  1) 


where  the  parameters  in  a  =  (xo,vx,ax,t/o,  vy,ay)  identify  a  particular  target  whose  track  is  defined  by 
p(t,a)  =  (x(t,a),y(t,a)). 

At  a  discrete  set  of  scan  limes  (tic}k~i  (  h  <  h  <  ...  <  f.v),  a  radar  located  at  the  origin  in  this 
Cartesian  space  observes  error  contaminated  ranges  and  angles  of  the  targets  in  the  observation  space  which 
is  a  circle  with  radius  R  centered  at  the  origin.  Some  observations  are  spurious  and  some  observations  of 
true  targets  are  missed.  At  time  tk,  the  radar  is  assumed  to  return  the  set  of  observations  {*,*  where 
Mk  is  the  number  of  observations  and  z*k  =  (rt*  ,  0*  ).  To  every  scan,  a  dummy  observation  Zq  is  added  to 
represent  missed  detections.  Each  observation  (r*  ,  0*  )  arising  from  an  existing  target  is  related  to  the  true 
observable  H(p(tk,ac ))  by 

(c.2, 
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where  ek  and  ek  are  independent  zero-mean  Gaussian  random  variables  with  standard  deviations  a*  = 
ar(£*;)  and  cx£  —  ag(lk),  respectively.  The  measurement  error  covariance  matrix  is  given  by  5vr«(£)  = 
diag(a2(t), (More  generality  is  obtained  by  allowing  aT(l)  and  a6{t)  to  vary  with  the  spatial  position 
of  the  measurement.)  The  true  observable  H(p(t,a)}  is  related  to  the  track  p{t,a)  by 


ll(p{t,a)) 


y/x(t,a)2  +y(t,a)2\ 

arclan[|?TS)]  /' 


(C.3) 


(If  the  definition  of  arctan  is  based  on  the  principal  angle,  then  the  appropriate  shifts  in  must  be  made.) 
If  the  observation  corresponds  to  a  false  alarm  or  new  target,  then 


*k 

0k 


if  the  observation  is  spurious, 

if  the  observation  arises  from  a  new  source, 


(C.4) 


where  the  random  sequences  wk  and  uk  have  some  assumed  densities  pk  and  pk,  respectively.  A  common 
assumption  is  that  (u and  (uk,  uk)  are  uniformly  distributed  over  the  observation  space  so  that 


p*(wk)  = 


if  0  <  wk  <  R  and  —  7r  <  <  n 

0  otherwise 


(C.  5) 


A  similar  expression  would  hold  for  pk.  The  number  of  false  alarms  and  new  targets  are  assumed  to  be 
generated  at  each  time  interval  [£fc_ ) .  £*j  according  to  a  Poisson  distributions  with  expected  numbers  A*  and 
A*,  respectively. 

A  set  of  observations  Zu..lN  —  {z,1 ,  •  •  •  ,z?N },  containing  one  observation  from  each  scan,  will  be  called 
a  track  of  observations.  Note  that  some  observations  in  a  track  of  observations  might  be  dummy,  representing 
a  missed  detection.  Tracks  of  observations  f?o  oifco  0  with  a  single  nonzero  index  (i.e.,  a  single  non-dummy 
observation)  will  be  taken  to  represent  false  reports  (clutter). 

To  determine  the  most  probable  partition  of  the  observations  into  tracks  and  false  reports,  each  track 
of  observations  Zly  lN  must  have  a  likelihood  Ltl  ,N  associated  with  it.  First,  define  Lo-  o  —  1  and 
L 0  •  Oiko.  0  =  1-  Next,  if  Z*,.  has  at  least  two  nonzero  indices,  the  following  likelihood  expression  can  be 
derived  [37]: 


Ao‘*  fr(i 


p})p£p 


a yf(zkk\z,r  lN) 


fA^(2*jZ„...,J 

v) 


*} 


(1-AoiJ 


(C.6) 


where 


Pk 

r<t> 


Pk,  if  track  Z,,. ..tfl  terminates  at  scan  k\ 

(1  -  Pk)(  1  -  P$ ),  if  track  ZM ...,w  has  a  missed  detection  on  scan  k; 
1,  otherwise, 


and  the  indicator  functions  vk,  6k,  and  are  defined  by 


if  2*  is  a  new  target; 
otherwise; 

if  zk  belongs  to  an  existing  track, 
otherwise; 

'  =  j; 

otherwise. 


(C.7a) 


(C.7b) 
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Also,  P*  and  denote  the  payabilities  of  termination  and  and  detection  on  scan  k.  The  likelihood 
functions  pf(zXk \ZXl ..  lN  ),  p*(2,*  |E,,  .*  ),  and  p£(2*  | Z»,  t,v  J  are  those  ,,f  the  error  in  'he  observations  of  the 
target,  false  report  model,  and  new  source  model,  respectively. 

If  the  target  dynamics  (C.l)  are  known,  the  likelihod  expression  in  (C.6)  can  be  computed  using  the 
densities 

Pt(Zy =  - 7T~~T  -  77"-  exp  {  -  l-  [2*  -  //fc(p(£fc,a))]7  [2*  -  Uk{p{tk,  Or))]  }  (C.  8) 

2nar(tk)ag{tk)  2 

In  the  context  of  formulating  the  data  association  problem,  however,  the  track  parameters  are  unknown  and 
must  be  estimated  from  the  observations  {z^ 2*}.  The  parameters  a  are  replaced  by  the  maximum 
likelihood  estimate 

a  =  Arg  Max  L„.  ,J.p(-,a))  =  Arg  Min  c,,.  ,„(?(•, a)),  (C.9) 

which  can  be  equivalently  characterized  as  the  solution  to  the  nonlinear  least  squares  problem 
d  =  Arg  Min  c,,  (p(-,  a)) 


=  Arg  Min'Y'Xl  -  A 0<J  (2*  -  //*(p(£*.  a))]T  ErV(£fc)  (*,*  -  Hk(p(tk,a r))] 


(C.10) 


*:=  1 


III.C.2.  The  Data  Association  Problems.  This  section  will  address  the  formulation  of  the  data 
association  problem  for  two  cases,  in  track  initiation,  no  tracks  are  assumed  known.  In  track  maintenance, 
some  tracks  may  be  known  from  prior  information.  These  tracks  must  be  extended  using  the  newly  arriving 
information,  while  keeping  in  mind  that  new  tracks  might  also  be  initiating.  We  first  address  track  initiation. 
For  every  track  of  observations  Ziv  iN,  define  a  0-1  variable  2,,  ,.iN  via 

_  _/l  observations  { z], , ....  2*  }  were  generated  by  the  same  target  frilaf 

"  \  0  otherwise 

and  a  score  by 

c,,  ...,*  =  —  -is,  (CM  16) 

where  Lix...iN  was  defined  in  (C.6).  The  requirement  that  a  single  non-dummy  observation  2*  (1  <  k  <  N, 
1  <  »*;  <  Mk)  be  either  a  false  report  or  assigned  to  exactly  one  track  can  now  be  expressed  as 

M\  Mk- \  Mk-i  Mn 


E-  E  E  V,.. 

ti  =0  i*-i=0ik.i-0  t/v-0 


■It-  'N 


1. 


(C.12) 


The  data  association  problem  of  partitioning  the  observations  into  tracks  and  false  reports  can  now  be 
posed  as  the  following  multidimensional  assignment  problem 


M, 

mn 

Minimize 

E- 

]  CM  \s  zu  -is 

-0 

t/v  -0 

Mj 

Ms 

Subject  To  : 

E 

^  ^  Z\\  -is  1,-. 

0  -n 

IN  0 

Mi 

Mtr  t  Mtt*  1  Ms 

E 

E  E  -E 

u  0 

1  k  -  1  —  0  t  *  *- 1  •-  0  t/V-0 

for 

i*  =  1, . . . ,  Mk  and  k  =  2 

M 1 

Mn  -  , 

E- 

^  z»i  -  -\s 

m  -o 

*/v  -  i  -0 

z>i  »« 

€  {0,  1 )  for  all  »i, .. .  ,iy, 

(0.13) 
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Next,  track  maintenance  using  a  sliding  window  is  developed.  Suppose  ihat  the  observations  on  I’ 
previous  scans  (of  observations)  have  been  partitioned  into  tracks  and  false*  alarms  and  that  K  new  scans  of 
observations  are  to  be  added.  One  approach  to  solving  the  resulting  data  association  problem  is  formulate 
the  problem  as  a  track  initiation  problem  with  P-rK  scans.  This  is  the  previously  mentioned  batch  approach 
The  deferred  logic  approach  adopted  here  is  to  treat  the  track  extension  problem  within  the  framework  of  a 
window  sliding  over  the  observation  sets.  First  assume  that  the  scans  of  observations  arc  partitioned  into 
three  components:  D  discarded  scans  of  observations,  II  retained  scans  of  observations  from  the  P  previously 
processed  scans,  and  K  new  scans  of  observations.  Thus  the  number  of  scans  in  the  sliding  window  is 
jV  —  R  +  K  w*hile  the  number  of  discarded  scans  is  D  -  P  -  R. 

Let  A/o  denote  the  number  of  confirmed  tracks  previously  constructed  from  the  discarded  and  retained 
regions  that  are  present  at  the  start  of  the  tracking  window.  (These  Mr,  tracks  may  be  obtained  as  the 
solution  of  a  previous  problem  assignment  problem,  may  be  union  of  the  best  K  such  solutions,  or  may 
be  all  of  the  feasible  tracks.  However,  tracks  terminated  in  the  discarded  region  are  generally  not  included 
in  Mo.)  Suppose  the  iff1  such  track  is  denoted  by  T,0  for  io  =  l,...,Mo.  For  io  >  0,  the  (.V  —  l)-tuple 
{r,01  z,1,  .... ,  zff  }  will  denote  a  track  T,0  plus  a  set  of  observations  or  measurements  {s,1, , . . . ,  z?N  }  ,  actual 
or  dummy,  that  are  feasible  with  the  track  Tl0.  The  (  V  +  l)-tuple  { Tq , zf, . . ■ ,  zff }  will  denote  a  track  that 
initiates  in  the  sliding  window 

Analogous  to  the  track  initiation  case,  one  can  define  the  zero-one  variable 

1  if  {Tl0,< . **  }  is  assigned  as  a  unit,  (C.14a) 

0  otherwise. 

and  the  corresponding  cost  for  the  assignment  of  the  sequence  {TJ#,  ,  zXN  }  to  a  track  by 

c»o>i ■■■is  —  —  Iti  Lt0 L,,  ...lN .  (C.  146) 

Here  Lt,q  is  the  composite  likelihood  from  the  discarded  scans  just  prior  to  the  first  scan  in  the  window 
for  io  >  0,  Lt„  =  1,  and  LXy..tN  is  defined  as  in  (C.6)  for  the  N- scan  window.  (LTa  =  1  is  used  for  any 
tracks  that  initiate  in  the  sliding  window.)  The  data  association  problem  for  track  maintenance  can  thus  be 
formulated  as 

Mo  M/v 

Minimize  £  £  <V-.v  *»«••■»/* 

10  =0  ijv=0 

Mi  Mn 

Subj.  To  Y'...  Y  =  1.  M0, 

1 1  =0  t,v  =0 

Mo  M2  M/V 

^  ^  ^  ^  ■  ‘  ^  ^  ^to* i  1 

lo  -0  —  0  »,v  —0 

Mt  Mir  l  Mlc^-J  M/v 

E  -  Z  E  ■  E-o  -,v  =  i. 

*o-0  i/v-0 

for  ik  —  1 , . . . ,  M.v  and  k  —  2, . . . ,  .V  -  1 , 

Mo  M,v  .  t 

Z  "  E  z,° -■•n  ~  ~  . ^*v’ 

*o  -0  */v  -  i  -  0 

€  {0,1}  for  all  io,  ■  •  • ,  i,v- 

Note  that  the  association  problem  involving  N  scans  of  observations  is  an  ,'V-dimensional  assignment  problem 
for  track  initiation  and  an  (/V  +  1  (-dimensional  one  for  track  maintenance. 


(C.15) 
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III.C.3.  Numerical  Simulations.  The  performance  oi  the  tracking  algorithms  presented  in  this 
work  depends  on  many  factors,  including  target  density  (i.e.,  the  number  of  targets  per  unit  space),  the 
space  size,  the  measurement  error  covariances  and  probabilities  of  the  detection  and  false  report  rates.  This 
section  presents  results  of  two  studies  investigating  the  impact  of  changing  the  sensor  error  characteristics 
and  window  sizes  on  the  performance  of  these  algorithms.  Results  of  more  comprehensive  parametric  studies 
will  be  presented  in  a  future  work. 

The  tracking  problems  considered  in  both  studies  had  the  following  characteristics:  the  observation 
space  is  circular  with  a  radius  of  20  miles  and  sensor  in  the  center  and  10  targets  that  were  initiated  before 
the  first  scan  and  never  terminated  (targets  were  generated  so  as  to  r.  er  leave  the  space).  The  initial  target 
speeds  were  between  200  and  900  miles  per  hour  and  the  target  accelerati  o  was  not  more  than  0.0034  miles 
per  second  squared  (or  14,06-1  mill's  per  hour  squared).  The  scan  limes  were  every  10  seconds,  and  fifteen 
scans  of  observations  were  used.  .Note  that  the  total  gain  .n  velocity  due  to  acceleration  was  thus  not  more 
than  1500  mph.  The  sensor  returned  measurements  in  polar  coordinates,  as  described  in  Section  I1I.C.2.  All 
the  observations  in  each  scan  were  synchronous.  The  range  :rror  was  relative,  the  angle  error  was  absolute, 
and  neither  varied  with  time.  Both  the  missed  detections  and  false  reports  were  allowed  by  the  sensor. 

The  problems  were  generated  randomly  To  obtain  a  meaningful  sample  for  the  comparisons,  a  set  of 
100  problems  for  each  set  of  parameters  was  generated,  differing  only  in  the  random  number  generator  seed. 
All  results  presented  in  this  section  are  thus  averages  over  100  prot'ems. 

The  first  study  investigated  the  impact  of  changing  the  standard  deviations  of  the  measurement  errors. 
All  the  problems  in  the  first  study  had  the  probability  of  detection  of  0.95,  and  the  probability  of  false 
reports  (i.e.,  the  probability  that  an  arbitrary  observation  was  false  report)  of  0.05.  A  six  scan  window  was 
used  for  tracking.  Thus  a  six  dimensional  assignment  problem  governed  the  data  association  problem  for 
track  initiation  and  a  seven  dimensional  one  for  track  maintenance.  The  data  association  problems  arising 
from  the  tracking  problems  in  this  study  had  between  350  and  5000  variables  and  were  on  the  average  solved 
in  less  than  half  a  second  for  the  biggest  problems,  using  an  IBM  RS/6000-550  workstation.  Table  1  shows 
the  quality  of  tne  computed  tracking  problem  solution  for  each  of  nine  combinations  of  the  range  and  angle 
errors.  The  solution  quality  was  checked  after  each  scan  in  the  simulations.  Each  scan  thus  corresponds  to 
a  column  in  the  table.  The  column  for  6  scans  refers  to  track  initiation.  All  other  columns  are  results  for 
track  maintenance.  The  angle  error  standard  deviation  og  is  expressed  in  degrees. 
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Table  l  :  Solution  quality  for  varying  measurement  error 

Table  l  shows  that  the  quality  of  the  solution  increases  as  more  information  becomes  available  and  the 
track  estimates  become  better.  Six  scan  window  seems  adequate  for  tracking  problems  of  this  level  of  noise, 
especially  since  our  studies  of  individual  problems  indicate  that  the  quality  criteria  described  in  Section  4 


are  a  little  too  stringent  in  I  he  sense  tha’  tracks  proclaimed  "missed"  often  lie  just  outside  of  the  region 
defined  by  the  quality  criteria. 

For  the  second  study,  the  measurement  error  standard  deviations  were  kept  constant  at  ar  —  0.01  and 
=  0.5°,  while  the  probability  of  detection  was  either  0.!)5  or  0.7  and  the  probability  that  an  arbitrary 
report  was  fa!.>e  was  either  0.05  or  0.3.  For  each  of  the  four  combinations  of  these  two  parameters,  tracking 
problems  were  solved  using  I,  5,  fi,  and  7  scan  moving  windows.  Assignment  problem  sizes  ranged  from  30 
to  2500  variables  and  longest  solution  times  were  0  5-J  seconds  for  largest  windows  and  highest  noise  levels 
; most  of  the  cases  had  the  solution  times  of  less  than  0.3  seconds).  Table  2  shows  the  solution  quality  for  all 
16  combinations  of  probability  of  detection  (I'D),  probability  of  false  reports  f  1’FR)  and  window  size  ( Win  j 
inapplicable  entries  in  Table  2  are  denoted  by  **-” .  As  in  Table  1,  the  results  are  presented  on  scan  bv  scan 
basis,  with  the  first  number  in  each  rew  indicating  track  initiation,  while  the  remaining  numbers  refer  to 
track  maintenance. 
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Fable  2  Solution  quality  for  different  PD,  PFR  and  window  sizes 


Table  2  shows  that  the  quality  of  the  solution  again  increases  as  more  information  becomes  available. 
However,  in  problems  with  the  low  probability  of  detection,  even  after  15  scans  of  information,  the  algorithm 
that  uses  only  four  scans  in  the  moving  window  fails  to  identify  over  half  the  tracks.  As  the  window  increases, 
so  does  the  quality,  since  more  information  is  available  to  the  tracking  algorithm.  The  capability  to  vary  the 
window  sizes  is  thus  crucial  if  the  algorithm  is  to  handle  problems  with  different  noise  levels  successfully. 
Even  though  the  model  initiates  all  the  targets  before  the  lirst  scan  is  made,  the  results  in  Table  2  show  that 
track  initiation  must  fie  allowed  on  later  scans,  as  outlined  in  Section  3.  Especially  in  problems  with  low 
probability  of  detection,  some  targets  will  riot  be  indentilied  in  the  first  few  windows,  simply  because  they 
have  not.  been  detected  enough  times,  arid  are  thus  initialed  later. 

Comparing  Tables  1  and  2,  it  is  obvious  that  varying  the  probability  of  detection  has  more  impact  on 
the  solution  quality  than  varying  the  false  report  rate  or  the  measurement  rate.  This  is  not  really  surprising 
since  lowering  the  probability  of  detection  actually  removes  information  from  the  problem,  while  increasing 
the  false  report  rate  and  increasing  the  measurement  errors  just  adds  noise.  However,  using  the  appropriate 
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window  sizes,  the  algorithms  presented  in  this  work  construct  high  quality  solutions  even  for  very  noisy 
tracking  problems. 

Finally,  it  should  be  pointed  out  that  even  the  limited  parametric  studies  presented  m  this  section 
involved  solving  tracking  problems  with  widely  differing  amounts  and  types  of  noise  The  only  adjustment 
made  to  the  algorithms  was  the  size  of  the  sliding  window  Tins  indicates  that  the  algorithms  for  solving 
tracking  problems  using  multiscan  sliding  windows,  maximum  likelihood  estimation,  and  l.agrangian  relax¬ 
ation  for  data  association  problems,  are  quite  robust  and  thus  likely  to  be  effective  for  a  wide  range  of 
tracking  problems. 

III.  D.  Algorithm  Overview.  A  primary  objective  of  this  work  has  been  the  development  of  algorithms 
for  the  fast  construction  of  high  quality  :  near-optimal )  suboptirnal  solutions  of  the  following  multidimensional 
assignment  problem  fC.  13).  These  assignment  problems,  as  developed  in  Section  III .C,  possess  the  following 
important  characteristics:  the  problem  is  large  scale:  the  objective  function  is  noisy  due  to  plant  noise,  errors 
in  the  sensor  measurements,  and  modeling  uncertainty:  it  is  NT- hard 12  but  must  be  solved  real-time.  Gating 
and  clustering  techniques  are  generally  used  to  reduce  the  size  and  complexity  of  the  problem,  thereby  making 
the  problems  sparse.  We  argue  that  the  problem  should  be  solved  to  the  noise  level  and  not  to  optimality, 
since  the  objective  is  to  use  this  assignment  problem  as  a  vehicle  to  identify  objects  in  sensor  fusion  and 
estimate  tracks  in  tracking  The  NP-hardness  and  real-time  needs  rule  out  conventional  techniques  such  as 
branch  and  bound  or  explicit  enumeration.  The  sparsity  of  the  problem  raises  the  issue  of  whether  or  not 
(D.l)  has  a  feasible  solution.  To  resolve  this  we  assume  that  all  zero-one  variables  with  exactly  one  nonzero 
index  are  free  to  be  assigned  and  the  corresponding  cost  coefficients  are  defined.  The  zero-one  variable 
z0  o  and  the  corresponding  term  cq  ,>*0  o  in  the  objective  function  are  present  for  notational  convenience. 
Finally,  other  problems  of  interest  include  the  situation  in  which  the  “  =  1"  in  the  constraint  (D.l)  is  changed 
to  “  <,  =,  or  >  n*  for  some  nonnegative  in.  ^er  n*  However,  we  shall  not  address  these  problems. 

The  algorithm  development  in  this  work  is  based  on  Lagrangian  relaxation,  which  originally  gained 
prominence  as  a  method  for  efficiently  obtaining  tight  bounds  for  a  branch  and  bound  algorithm  in  Held 
and  Karp's  highly  successful  work  on  the  traveling  salesman  problem.  Overviews  of  this  methodology  can 
be  found  in  the  works  of  Geoffrion.  Fisher,  Shapiro,  the  book  by  Nemhauser  and  Wolsey,  and  the  references 
therein.  The  particular  algorithm  developed  in  this  work  is  motivated  by  that  of  Frieze  and  Yadegar  for 
three  dimensional  assignment  problems;  however,  an  overview  of  the  algorithms  developed  in  this  work  is 
perhaps  more  easily  described  interms  of  a  prototype  algorithm  for  a  general  integer  programming  problem. 

Consider  the  integer  programming  problem 

Minimize  c  z  =  V(z) 


Subject  To  .-fz  >  b 
Bz  >  d 


(D.  2) 


z,  is  an  integer  for  i  e  l  ■ 


where  the  partitioning  of  the  constraints  is  natural  in  some  sense.  The  l.agrangian  relaxation  of  (D .2)  relative 
to  the  constraints  Bz  >  d  is  defined  to  be 


<f*(u)  —  Minimize  o(z,u)  =  {crz  —  u1  (Bz  —  d)} 
Subject  To  Az  >  b 


(D.  3) 


z,  is  an  integer  for  i  €  / 
u  >  0 
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where  u  >  0  is  interpreted  componentwise.  It  the  constraint  set  Bz  s  d  is  replaced  bv  Hz  —  d.  the 
nonnegativity  constraint  on  u  is  removed.  C  =  c1  z  -  u!  ( Bz  —  d\  is  the  Lagrangian  relative  to  the  constraints 
Bz  >  d,  and  hence  the  name  Lagmngtan  relaxation.  Next,  if  z  is  an  optimal  solution  to  il),2:.  problems 
(D.2)  and  (D.3)  imply 

for  ali  a  >  0  \DA) 

Given  a  specific  multiplier  u,  let  zr{u)  and  rriul  denote  suboptirnal  and  optimal  solutions  of  the  relaxed 
problem  (D.3),  respectively.  Generally,  zr(n;  is  not  feasible  for  the  relaxed  constraint  Bz  >  d:  however,  if 
zr(u)  is  feasible,  then  it  is  then  also  optimal  for  (D.2).  Thus  one  must  develop  a  recovery  procedure  for 
constructing  a  feasible  solution  of  (D.2)  from  either  of  zT(u)  or  ;r(u).  There  are  several  reasons  why  the 
resulting  feasible  solution  zj  might  be  a  good  solution  of  (D.2).  First,  if  the  multiplier  u  =  a  >  0  is  chosen 
as  the  maximizer  of  the  problem 

Maximize  (4>(u)  :  u  >  0}  (O.i) 

and  the  duality  gap  [<t>(u),  V^i)]  is  small,  then  the  recovered  feasible  solution  zj  of  (D.2)  from  the  solution 
iT(ti)  of  (D.3)  may  be  close  to  z.  (The  experience  of  many  researchers  is  that  this  duality  gap  tends  to  be 
much  smaller  for  equality  constrained  problems  than  for  corresponding  inequality  constrained  problems :iJ.) 
Secondly,  the  term  -ur(Bz  -  d)  in  (D.3)  acts  like  a  penalty  for  violating  the  constraint  Bz  -  d  >  0,  thereby 
forcing  zr(u)  closer  to  the  optimal  solution  z  of  (D.2).  Finally,  the  recovery  procedure  should  be  designed  to 
minimize  any  remaining  flexibility  in  the  objective  function  in  (D.l).  Thus  given  this  rationale,  the  following 
prototype  algorithm  abstracts  some  of  the  ideas  of  the  work  on  three  dimensional  assignment  problems  by 
Frieze  and  Yadegar: 

Prototype  Algorithm.  Construct  a  sequence  of  multipliers  {u*}£°:0  in  the  course  of  maximizing  <F(u) 
defined  in  (D.3)  and  a  corresponding  sequence  of  feasible  solutions  {z^}^  0  of  (D.2)  as  follows : 

A.  Choose  an  initial  approximation  n0. 

B.  Given  u*,  determine  a  new  multiplier  u/t-i  from  a  step  in  the  maximization  problem  (D.5),  so  that 
$(ufc)  <  'I'ko). 

C  Given  ttu and  a  solution  zr(uic .  i )  of  (D.3),  recover  a  feasible  solution  r* . 1  of  the  integer  programming 
problem  (D.2). 

In  the  absence  of  any  a  prior  knowledge  of  the  initial  multiplier  no,  a  good  neutral  choice  in  Part  A  is 
u0  =  0.  Part  B  of  this  algorithm  is  the  nonsmooth  optimization  phase  and  one  of  the  most  widely  used 
methods  for  non-smooth  optimization  is  the  subgradient  algorithm,  which  is  the  nonsmooth  analog  to  the 
steepest  ascent  method.  Analogous  to  conjugate  gradient  methods  for  smooth  optimization  is  the  class  called 
“bundle  methods”.  This  includes  the  space  dilation  method  of  Shor.  the  “bundle-trust  region”  method  due 
to  Schramm  and  Zowe,  and  the  conjugate  subgradient  method  of  Wolfe.  Wolfe’s  algorithm  is  used  in  this 
work.  The  recovery  procedure  is  part  C  of  this  algorithm  and  can  vary  considerably  with  the  problem.  Note 
that  <f>(ufc)  <  <F(u)  <  V(z)  <  V{zk),  so  that  we  have  a  bound  on  the  optimal  solution.  With  an  estimate  of 
the  noise  level  in  the  problem,  we  can  then  use  these  bounds  as  a  stopping  criteria.  The  explicit  developed 
for  the  multidimensional  assignment  problems  are  presented  in  [2d]  and  forthcoming  work. 
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IV.  Infinite  Dimensional  Optimization  Problems  and  Numerical  Control 

The  classes  of  optimal  control  problems  currently  under  investigation  are*  subclasses  of  the  problem 

Minimize  Jjx,  u}  :=  ^Uo,  i * .  x(tt ))  -r  j  /,»(£,  x.  u  ult 

Subject  To  i  —  f{t,x[t),u(L)) 

/i(£o,fi.x(fo).x(ti))  =  0 
h(t,  x{t),u(l))  ~  0  a.(*.  [£(>,£ ij 
<j(t,  x(t),u[L))  <  0  a.e.  [to.£i] 
u(t)  €  V.  a.e.  (£0l  ‘ ,] 

jT,(£o.  1 1, x(£o). z(£i ), x,  u)  <  0  for  t  =  1 , . . . ,  k 

(X,u)€  ^•“(Mtl.nt")  X  L°°{[t0,tx}Mm) 

where  x  is  an  n- vector,  u  is  an  m- vector,  11  is  a  boundary  operator,  Q  is  a  closed  convex  set,  Tj  is  a 
functional  and  Wl,p([to,  ti  ],  HU1)  is  the  usual  Sobolev  space  which  can  be  characterized  via  the  Sobolev 
imbedding  theorem  as  consisting  of  those  absolutely  continuous  vector  functions  with  the  first  derivative  in 
^([io,  1 1],  IR").  The  functions  ,r,  /0,  /,  B ,  /i,  Tx,  and  g  are  assumed  to  be  at  least  C2  with  respect  to  their 
arguments. 

Our  interest  in  this  problem  is  two  fold.  First,  working  with  former  PhD  student  B.  Yang,  Professor 
VV.  W.  Hager  of  the  University  of  Florida  and  Professor  Asen  Dontchev  of  Bulgiarian  Academy  of  Sciences 
and  Math  Reviews,  we  have  investigated  the  convergence  of  various  numerical  methods  (Newton’s,  penalty, 
augmented  Lagrangian,  interior  point  methods)  in  the  appropriate  infinite  dimensional  spaces.  This  work 
has  evolved  as  follows:  Poore,  Yang,  and  Hager  [42]  have  investigated  convergence  of  penalty,  multiplier,  and 
Newton  methods  for  a  subclass  of  the  above  problems  without  set  constraints  and  inequality  constraints  on 
the  controllers.  A  more  theoretical  analysis  of  the  above  problem,  again  without  set  constraints  or  inequality 
pointwise  constraints  on  the  state  variables  and  controller  was  developed  in  the  PhD  thesis  of  B.  Yang 
[50].  This  latter  work  has  been  considerably  generalized  to  include  these  pointwise  equality  and  inequality 
constraints  in  the  recent  work  of  Dontchev,  Hager,  Poore,  and  Yang  [11].  The  approach  was  to  first  derive 
sufficient  optimality  conditions  for  an  infinite  dimensional  optimization  problem  in  a  setting  that  is  applicable 
to  optimal  control  problems  with  endpoint  constraints  and  with  equality  and  inequality  constraints  on  the 
controls.  Under  the  hypotheses  of  the  sufficient  optimality  theorem  we  show  that  the  solution  to  an  optimal 
control  problem  possesses  a  Lipschitz  stability  property  relative  to  problem  perturbations.  As  an  application 
of  this  stability  result,  we  establish  convergence  results  for  the  successive  quadratic  programming  algorithm 
and  for  penalty  and  multiplier  approximations  applied  to  optimal  control  problems. 

The  second  area  of  research  interest  is  the  parametric  problem  associated  with  the  above  optimal  control 
problem.  The  interaction  of  multiple  and  bifurcating  states  in  the  absence  of  controls,  periodic  phenomena, 
chaotic  behavior,  and  bifurcating  controls  arising  from  the  dynamical  systems  and  holonomic  constraints  is 
open  to  investigation.  Given  a  certain  phenomena  arising  from  a  dynamical  system,  the  problem  may  be  to 
control  this  phenomena,  to  determine  multiple  solutions,  or  to  investigate  the  dependence  of  a  solution  on 
the  system  parameters  over  a  wide  range,  i.e.,  global  sensitivity.  (The  latter  is  also  important  in  adaptive 
control.)  The  development  and  use  of  theoretical  and  numerical  bifurcation  and  continuation  methods  in 
dynamical  systems  and  nonlinear  equations  has  been  spectacularly  successful  in  analyzing  and  understanding 
the  phenomena  represented  by  these  systems,  but  we  know  of  no  systematic  treatment  or  works  on  the 
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constrained,  nonlinear  parametric  control  problem  paralleling  that  found  in  dynamical  systems  Thus  a  long 
term  goal  of  this  research  program  will  be  the  investigation  of  parametric  problems  in  nonlinear  control 
systems  including  but  not  limited  to  the  nonlinear  optimal  control  problem. 

V.  Parametric  Nonlinear  Programming  and  Control 

This  section  describes  the  parametric  optimization  problem,  some  of  the  accomplishments  over  the  last 
two  years  and  an  application  in  design  optimization. 

V.A  Problem  Statement.  The  use  of  bifurcation  and  singularity  theory  in  the  investigation  of  para¬ 
metric  problems  in  optimization  and  control  represents  a  potential  for  a  real  extension  of  our  understanding 
of  basic  phenomena,  global  sensitivity,  robustness,  and  multiplicity  of  solutions  in  both  finite  and  infinite 
dimensional  optimization  problems,  in  much  the  same  wav  that  these  theoretical  and  numerical  techniques 
have  helped  our  understanding  of  dynamical  systems  and  nonlinear  equations.  Thus  the  objective  in  this 
program  has  been  the  development  of  the  analytical  and  numerical  techniques  to  map  out  regions  of  qual¬ 
itatively  different  behavior  and  to  locate  the  “stability”  boundaries  of  these  regions  in  parameter  space. 
The  latter  is  important  because  drastic  changes  in  the  optimum  occur  at  singularities,  which  define  these 
“stability”  boundaries.  Our  initial  work  (28, 35, <17, -18)  has  been  the  classification  and  analysis  of  singular 
points  in  the  nonlinear  parametric  programming  problem.  Numerical  techniques  for  predictor-corrector  con¬ 
tinuation  techniques  have  been  developed  using  a  nonstandard  variable  order  Adams-Bashforth  predictors 
with  an  adaptive  error-step  size  control  [19j.  This  software,  which  is  available  through  Netlib,  is  particularly 
efficient  and  robust  for  parametric  problems.  Numerical  continuation  and  bifurcation  techniques  are  being 
developed  and  tailored  to  the  finite  dimensional  constrained  optimization  problem  in  support  of  future  work 
on  large  scale  control  systems  design  optimization  (20, 23, 2*1).  In  the  remainder  of  this  section,  the  paramet¬ 
ric  nonlinear  programming  problem  is  defined  and  some  of  the  accomplishments  during  1991  and  1992  are 
presented. 

Mathematically,  the  parametric  nonlinear  programming  problem  is  that  of  determining  the  behavior  of 
solution(s)  as  a  parameter  or  vector  of  parameters  a  €  IRr  varies  over  a  region  of  interest  for  the  problem 

Minimize  {f(x,  a)  I  c,(x,a)  =  0  for  i  £  E 

{A.l) 

c,(x,  a)  <  0  for  i  £  /} 

where  E  —  {1, . . .  ,p}  and  /  =  {p  4-  1, . . .  ,p  +  q)  represent  the  index  sets  for  the  equality  and  inequality 
constraints,  respectively,  and  where  /  :  Uln~r  — *  IR,  C£  :  IRn"r  —  1R9  and  c/  :  Hln^r  — ►  W  are  assumed  to 
be  at  least  twice  continuously  difTercntiable.  Using  first-order  necessary  conditions  as  motivation,  one  can 
convert  the  characterization  of  a  solution  to  this  problem  as  a  solution  of  a  system  of  nonlinear  equations. 
At  a  regular  point  of  this  latter  system,  the  implicit  function  theorem  rigorously  justifies  the  computation  of 
the  derivatives  of  the  primal  and  dual  variables  with  respect  to  the  parameter  a.  These  derivatives  provide 
the  basis  for  local  sensitivity  analysis  as  presented  in  the  work  of  Fiacco  [12,13]  and  references  therein.  Thus 
all  the  “action"  is  at  the  singular  points  of  this  system  (A.l)  where  catastrophic  failure,  extreme  sensitivity, 
and  jumps  to  undesirable  operating  stales  can  occur.  We  have  investigated  these  singularities  in  several 
papers  [28,35,47,48]  via  bifurcation  and  singularity  theory.  (These  singular  points  are  characterized  by  a  loss 
of  strict  complementarity,  a  violation  of  the  linear  independence  constraint  qualification,  or  the  singularity 
of  the  Hessian  of  the  Lagrangian  on  the  tangent  space  to  the  active  constraints.) 

The  work  in  the  last  two  years  has  turned  to  the  development  of  numerical  continuation  and  bifurcation 
techniques  for  the  systematic  exploitation  of  these  methods  in  applied  problems.  We  now  give  a  synopsis  of 
the  algorithms  and  methods  that  can  be  found  in  a  scries  of  papers  by  Poore  and  Lundberg  [19-24]. 
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A  solution  of  the  parametric  programming  problem  (A.  1 )  is  a  solution  of  the  following  system  of  nonlinear 
equations 


F(x,  A,  u\ a)  = 


'vxC(x,\,w,a) 
A c(r;  q) 
v2  +  ArA  -  :32 


=  0 


(--1.2) 


where  £  =  C(x,  \,u)  =  e/(i)  +  A,c,(x)  is  the  Lagrangian  and  A  is  a  diagonal  matrix  with  Alt  =  1 

for  i  e  E  and  A„  =  A,  for  i  €  l .  Lot  A  =  />’ U  {t  €  /  :  e,(i,a)  —  0}  denote  the  active  set.  Note  that  any 
solution  of  this  system  at  which  u  >  0,  c,(x)  <  0,  and  Aj  >  0  for  all  i  e  /  satisfies  the  Fritz  John  first-order 
necessary  conditions.  This  system  also  employs  a  nonstandard  normalization  v2  +  A7  A  -  32  =0,  where 
/3b  is  a  fixed  positive  real  number.  The  standard  normalization  u  =  1  is  not  employed  since  it  requires  a 
constraint  qualification  for  iis  validity  and  the  violation  of  the  linear  independence  constraint  qualification 
is  a  singularity  in  the  above  system  (A. 2). 

Since  a  multiplier  corresponding  to  an  inactive  constraint  is  zero,  the  system  (A. 2)  can  be  reduced  in 
complexity  by  using  an  active  -set  strategy.  The  inactive  constraints,  i.e.,  those  c,  for  which  i  €  l  -  A.  are 
thus  removed,  yielding  the  active  set  system 


F(z,  a) 


X 

c(i,  a) 

B(\,u) 

=  0, 

where  z  = 

A 

V 

(A.3) 


m  =  n  +  \A\  +  1,  A  =  (Ai, . . . ,  Ap,  X,eA^t),  and  c  =  (cj, . . .  ,Cp,  c,€^n/),  £(z,a)  =  uf(x,  a)  +  YhtA  -X»c,(x,a) 
and  S(A,  i/)  =  v2  +  ArA  -  3 1-  Continuation  for  the  system  (A.3)  along  with  locating  the  zeros  in  one  or 
more  of  the  active,  inequality  multipliers  At,  i  €  Adi  or  an  inactive  constraint  c,  for  i  e  I  —  A  and  changing 
the  active  set  appropriately  is  then  equivalent  to  continuation  for  the  full  system  (A. 2). 


V.B  Status  of  the  Algorithms.  Numerical  algorithms  for  numerical  linear  algebra  in  the  continua¬ 
tion  procedure,  critical  point  type,  singularity  detection  and  classification,  and  branch  switching  have  been 
developed  in  three  papers  of  Poore  and  Lundberg  [20,23,24]  with  additional  work  in  progress.  We  now  give 
a  brief  overview  of  these  results. 

The  linear  systems  that  arise  in  the  continuation  steps  can  be  reduced  via  block  elimination  algorithms 
[6]  to  the  solution  of  several  linear  systems  of  the  form 


W 


Ax  1 

!  9 

_AA  ' 

M 

where 


A 

0 


(BA) 


the  matrix  H  is  the  Hessian  of  the  Lagrangian  or  some  approximation  to  it,  and  AT  =  Dxc(x,  a).  During  the 
continuation  steps  the  matrix  H  need  not  be  positive  definite  on  the  tangent  space  to  the  active  constraints; 
however,  both  null  and  range  space  methods  are  easily  modified  to  form  the  basis  for  the  linear  algebra 
steps(20j. 

The  classification  of  critical  point  type  is  based  on  [20] 


(B. 2a)  sign  i/, 

(B.2b)  signs  of  c,(x,  a)  for  i  €  l  -  A  and  A,  for  i  €  /  n  .4, 

(B. 2c)  signs  of  the  eigenvalues  of  V2Ct , 

where  V2Cr  denotes  the  restriction  of  the  Hessian  of  the  Lagrangian  to  the  tangent  space  of  the  active 
constraints  M(AT).  ft  is  only  the  latter  class  of  signs  (B.2c)  that  require  computation,  and  this  can  be 
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accomplished  by  computing  the  inertia  of  the  reduced  Hessian,  which  can  be  accomplished  by  either  null  or 
range  space  methods  [20J. 

Methods  for  detecting  singularities  due  to  the  loss  of  strict  complementarity,  loss  of  the  linear  indepen¬ 
dence  constraint  qualification,  and  singularity  of  the  Hessian  of  the  Lagrangian  on  the  tangent  space  to  the 
active  constraints  have  been  extensively  developed  in  the  work  of  Lundberg  and  Poore  [20] .  The  philosophy 
behind  singularity  detection  is  to  skip  over  them  during  the  continuation  procedure,  detect  their  presence, 
and  then  take  the  appropriate  action,  e  g.,  switch  branches,  change  orientation,  or  continue  along  the  current 
branch.  The  detection  and  classification  for  simple  bifurcations  and  folds  have  been  show  to  be  inexpensive 
‘by-products’  of  the  continuation  procedure  [20j;  however,  due  to  the  technically  detailed  classification  we 
omit  a  further  discussion  of  this  problem. 

V.C.  A  Model  Problem  from  Design  Optimization.  The  numerical  continuation  techniques  de¬ 
scribed  in  the  previous  sections  will  now  be  used  to  obtain  a  “global”  analysis  of  the  sensitivity,  stability, 
and  multiplicity  of  minima  for  a  parametric  nonlinear  programming  problem  arising  from  design  optimiza¬ 
tion.  The  problem,  w'hich  is  simple  yet  still  exhibits  the  basic  phenomena,  involves  the  design  of  a  two  bar 
planar  truss  with  semi-span  1,  unloaded  height  h,  and  load  p  as  indicated  in  Figure  1. 


Figure  1:  Loaded  Two  Bar  Truss 

Given  a  specific  unloaded  height  h  and  load  p,  the  deflection  d  is  a  minimizer  of  the  potential  energy 
E(d,h:p)  =  ~pd  -e  ^vT  -(-  h2  -  \J \  +  (h  -  d)2j  /vT  +  k2.  Rheinboldt  [44]  used  this  model  problem  to 
illustrate  continuation  methods  in  structural  analysis  and  has  given  a  rather  complete  solution  to  both  the 
static  and  parametric  problems.  Rao  and  Papalambros  [43]  posed  a  corresponding  optimal  design  problem 
as  that  of  choosing  the  height  h.  to  minimize  the  deflection  subject  to  0  <  h  <  1.5.  This  problem  is  posed 
mathematically  as 


Minimize  d 

(C.l)  Subject  To  V ^E{d,h\p)  =  0 

0  <  h  <  1.5 

In  addition  to  selecting  the  minimizer,  the  state  (d,  h)  must  also  be  selected  so  that  the  potential  energy 
E(d,  h;p)  is  minimized  with  respect  to  d.  The  corresponding  paiametric  problem  is  to  determine  the  solution 
and  its  properties  as  the  load  p  varies  over  all  physically  important  ranges.  The  numerical  methods  discussed 
in  the  work  of  Lundberg  and  Poore  and  briefly  discussed  in  the  previous  subsection  were  used  to  obtain  a 
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global  solution  of  this  problem  as  shown  in  Figure  2 


Figure  2:  Solutions  of  (A. 2)  for  Problem  (C.l) 

Here,  the  displacement  d  and  unloaded  height  h  as  p  varies  and  represents  a  projection  of  the  solutions 
of  (C.l)  into  ( h,d,p )  space.  Solid  and  dashed  lines  indicate  paths  of  local  minimizers  and  maximizers, 
respectively.  The  dashed  and  dotted  line  represents  a  feasible  singular  path,  and  lines  of  small  dots  represent 
infeasible  solutions.  The  solutions  to  the  optimization  problem  need  not  be  points  of  minimum  potential 
energy  E(d,  h;p),  which  is  not  minimized  on  the  segments  from  (b)  to  (c)  and  from  (d)  to  (c)  to  (e).  However, 
all  other  feasible  path  segments  do  correspond  to  physical  states  of  the  truss  where  the  potential  energy  is 
minimized. 

We  now  describe  these  singularities  and  the  connecting  path  segments,  beginning  with  those  which  occur 
along  the  solution  branch  where  the  constraint  h  <  1.5  is  active.  Loss  of  strict  complementarity  gives  rise  to 
the  bifurcation  points  (g),  (a),  and  (c),  whose  presence  was  indicated  by  a  change  in  sign  of  the  multiplier 
\a.  At  these  points  the  inequality  constraint  becomes  weakly  active  and  solution  paths  bifurcate  into  the 
region  0  <  h  <  1.5.  The  fold  points  (b)  and  (d)  (p  =  ±.37),  which  resulted  from  violation  of  the  linear 
independence  constraint  qualification,  were  detected  by  a  change  in  the  sign  of  v.  The  type  of  the  solution 
along  this  solution  branch  is  determined  by  the  sign  of  A0/iz,  which  changes  at  each  of  these  five  singular 
points.  This  results  in  the  alternating  segments  of  minimizers  and  maximizers  shown  in  Figure  2. 

The  solution  to  the  parametric  design  problem  can  now  be  described  for  p  >  0.  Given  a  small  but 
positive  load  p,  the  global  minimum  occurs  on  the  branch  of  minimizers  between  singular  points  (f)  and  (a). 
As  the  load  p  is  increased  from  zero,  the  height  h  increases  from  \/2  to  h  =  1.5  where  the  constraints  h  <  1.5 
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becomes  active.  As  the  load  p  is  increased  further,  the  deflection  d  continues  to  increase  along  the  path  from 
(a)  to  (b)  until  p  reaches  0.37  where  the  truss  ‘snaps  through’  and  there  is  no  minimum  beyond  p  =  0.37 
corresponding  to  a  height  h  near  1.5.  (The  only  way  to  maintain  an  optimum  locally  beyond  p  =  0.37  is  to 
increase  the  parameter  0  —  1.5  in  the  upper  bound  on  the  height  h.)  The  local  minimizer  corresponding 
to  h  —  0  becomes  the  global  minimizer  for  p  beyond  p  =  0.37,  Local  sensitivity  is  surely  present  at  points 
(a)  and  (b).  Note  that  the  path  of  minimizers  is  continuous  but  not  differentiable  at  (a).  (Near  such  points 
many  optimization  codes  exhibit  cycling.)  At  the  fold  point  (b),  the  path  of  minimizers  ceases  to  exist. 
Optimization  codes  would  have  difficulty  here  since  the  unnormalized  multipliers  will  be  large  and  go  to 
infinity  as  p  approaches  0.37.  The  conclusion  with  regard  to  the  design  of  the  truss  is  that  for  stability  the 
loads  must  be  less  than  p  —  0.37  and  that  sensitivity  occurs  near  the  singular  points  (a)  and  (b)  for  the 
reasons  stated.  Clearly,  the  ability  of  the  continuation  procedure  to  locate  such  singular  points  and  obtain 
such  a  global  analysis  is  a  major  strength  of  the  methodology. 
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2.  Multitarget  Tracking  and  Multidimensional  Assignment  Problem0,  SPIE  Meeting,  Orlando,  FL,  April, 
1991. 

3.  Parametric  Problems  in  Constrained  Optimization  and  Control,  University  of  Florida,  Gainesville,  FL, 
April,  1991. 

4.  The  Data  Association  Problem  in  Multitarget  Tracking  and  Multidimensional  Assignment  Problems, 
SDIO  Parallel  Programmers  Meeting,  Falcon  AFB,  Colorado  Springs,  April  30,  1991. 

5.  Multi-Target  Tracking  on  Parallel  Processors,  Space  Tech,  Fort  Collins,  Colorado,  May  20,  1991. 

6.  Multitarget/ Multisensor  Tracking  and  Multidimensional  Assignment  Problems,  US  Army  SDC,  Hun- 
stville,  Ala.,  June  10,  1991. 

7.  Multitarget  Tracking  and  Multidimensional  Assignment  Problems,  SDIO  Meeting  Hunstville,  Ala.,  June 
13,  1991 
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8.  Multitarget./ Multisensor  Tracking  and  Data  Association  Problems,  AFOSR/  Home  Lab  Workshop  on 
Multisensor  Fusion,  Rome,  NY,  June  IS,  1991. 

9.  Parametric  Problems  in  Constrained  Optimization  and  Control,  Conference  on  "Numerical  Optimization 
Methods  in  Differential  Equations  and  Control,”  NCSU,  Raleigh,  NC,  July,  1991. 

10.  The  Data  Association  Problem  in  Multitarget/Multisensor  Tracking,  IBM,  Boulder,  CO  and  Martin 
Marietta,  Denver.  CO,  August,  1991. 

11.  Numerical  Continuation  Methods  in  Parametric  Nonlinear  Programming,  International  Conference  on 
Parametric  Optimization  and  Related  Topics  Ill,  Giistrow,  Germany,  August,  1991 

12.  Some  New  Approaches  to  Multitarget/Multisensor  Tracking,  Army  Research  Office,  Research  Triangle 
Park,  November  15,  1991. 

13.  Numerical  Continuation  and  Bifurcation  Methods  in  Constrained  Optimization,  University  of  Florida, 
Jan.,  1992. 

14.  The  Data  Association  Problem  in  Multitarget/Multisensor  Tracking,  University  of  Florida,  Jan.,  1992 

15.  Multitarget/Multisensor  Tracking,  Rome  Labs,  Grilliss  AFB,  Feb  ,  1992 

16.  Multitarget/Multisensor  Tracking,  Federal  Sector  Division  of  IBM,  Feb.,  1992 

17.  Multitarget/Multisensor  Tracking,  DARPA,  Washington,  DC,  March  25,  1992 

18.  Data  association  for  track  initiation  and  extension  using  multiscan  windows,  1992  SPIE  Conference  on 
Small  Targets,  Orlando,  FL,  April,  1992. 

19.  Data  association  for  track  initiation  and  ex’ension  using  multiscan  windows,  Los  Alamos  Days  at  Boul¬ 
der,  CU-Boulder.  April,  1 992. 

20.  Data  association  for  multisensor  data  fusion  and  rnultitarget  tracking,  Grumman  Corporation,  Long 
Island,  NY,  July,  1992. 

21.  Data  association  for  multisensor  data  fusion  and  rnultitarget  tracking,  Rome  Labs,  Griffiss  AFB,  July, 
1992. 

22.  Data  association  for  multisensor  data  fusion  and  multitarget  tracking,  IBM  Corporation,  Owego,  NY, 
August,  1992. 

23.  Numerical  Continuation  and  Bifurcation  Techniques  for  Parametric  Nonlinear  Programming,  Fourth 
AI  AA/USAF/NASA/OAi  Symposium  on  Multidisciplinary  Analysis  and  Optimization,  Cleveland,  Ohio 
September  22,  1992. 


